#ifndef __radius_interpolation_cpp
#define __radius_interpolation_cpp

#include <iostream>
#include <fstream>
#include <iomanip>
#include "vect_3d.h"
#include "head.h"
//#include "class_ma.cpp"
// #include "class_ma_Complex.cpp"
//#include "class_mac16.cpp"
#include "svd_double.h"
// #include "svdlongdouble.cpp"
// #include "svdcmpld.cpp"
#include "gaulen.h"
// #include "cgmom.cpp"
#include <string>
#include <deque>
#include <ctime>
#include <limits>

using namespace std;

double lamda[299][300], kesai[299][300];

// lamda[][] is a 2D array used to store weight value of Gauss integral.

// kesai[][] is a 2D array used to store nodes value of Gauss integral.

// extern "C"

//{void __stdcall SVD(int *NRA,int *NCA,double *A,int *LDA,double *TOL,int *IRANK, double *GINVA,int *LDGINV);}

double* w, * x;

// double *P[100000];

// double *A[10000];

// double spx[100000];  //Scattered points.

// double spy[100000];

// double spz[100000];

dvector shift = dvector(0, 0, 0);

void generate_Gaussian_nodes() {

    x = new double[300];

    w = new double[300];

    int i, j;

    for (i = 0; i < 298; i++) {

        gaulen(-1., 1., x, w, i + 2);

        for (j = 0; j <= i + 1; j++) {

            lamda[i][j] = w[j + 1];

            kesai[i][j] = x[j + 1];
        };
    }

    delete []w;

    delete []x;
}

/*void input_scattered_points(char *nfile,int &num_n,float *spx,float *spy,float *spz)

{

   ifstream fin;

   char c;

   int i,num,num1;

   bool endf;







   ////////////////////////////////////////////////////////////////////////////////////////////////

   //2. create t_nlist.txt according to nlist.txt from ansys

   ////////////////////////////////////////////////////////////////////////////////////////////////





   fin.open(nfile);



   num=0;

   endf=false;

   for(i=0;i<50;i++)

   {

      fin>>c;

      if(c==NULL)

      {

         endf=true;

         break;

      }

   }



   while(endf==false)

   {

      for(i=0;i<7;i++)

      {

         fin>>c;

         if(c==NULL)

         {

            endf=true;

            break;

         }

      }



      if(endf==false)

      for(i=0;i<20;i++)

      {

         fin>>num1;

         if(num1==num)

         {

            endf=true;

            break;

         }

         else

         {

            num=num1;

            fin>>spx[num-1]>>spy[num-1]>>spz[num-1];

         }

      }

   }

   num_n=num;

   fin.close();



   ///////////////////////////////////////////////////////////////////////////////////////////////////////



}     */

void generate_new_interpolation_zeros(dvector* r, int N, double length) {

    double dl, dl_2;

    double* zeros, _pi = 3.1415927;

    int n, i, j, k;

    int section = 4;

    zeros = new double[2 * N];

    dl = length / (N - 1);

    dl_2 = dl / 2;

    for (k = 0; k < N; k++)

        // zeros[k]=length/2.*(1-cos((2.*k+1.)/(2.*N)*_pi));  //chebyshev nodes

        // zeros[k]=length/2.*(1-(cos((2.*k+1.)/(2.*N)*_pi))/(cos(_pi/(2.*N)))); //Li's chebyshev grid

        // zeros[k]=length/2.*(1-cos(k*_pi/(N-1)));

        // zeros[k]=length/2.*(1+kesai[N-2][k]);//gauss integral point

        // zeros[k]=length/2.*(1+tan((k+0.5)/(2.*N)*_pi-_pi/4.)); //tangent nodes

        zeros[k] = k * dl; // uniform

    n = 0;

    for (i = 0; i < N; i++)

        for (j = 0; j < N; j++)

            for (k = 0; k < N; k++) {

                // r[n]=dvector(i*dl,j*dl,k*dl);

                r[n] = dvector(zeros[i], zeros[j], zeros[k]);

                n++;
            }

    j = 0;

    for (k = 0; k < (N - 1); k++)

        // zeros[k]=length/2.*(1-cos((2.*k+1.)/(2.*(N-1))*_pi)); //chebyshev nodes

        // zeros[k]=length/2.*(1-(cos((2.*k+1.)/(2.*(N-1))*_pi))/(cos(_pi/(2.*(N-1)))));//Li's chebyshev grid

        // zeros[k]=length/2.*(1+kesai[N-3][k]); //gauss integral point

        // zeros[k]=length/2.*(1+tan((k+0.5)/(2.*(N-1))*_pi-_pi/4.));    //tangent nodes

    {

        zeros[j] = k * dl + dl / section; // uniform

        j++;

        zeros[j] = k * dl + dl * (section - 1) / section; // uniform

        j++;
    }

    for (i = 0; i < 2 * (N - 1); i++)

        for (j = 0; j < 2 * (N - 1); j++)

            for (k = 0; k < 2 * (N - 1); k++) {

                // r[n]=dvector(dl_2+i*dl,dl_2+j*dl,dl_2+k*dl);

                r[n] = dvector(zeros[i], zeros[j], zeros[k]);

                n++;
            }

    cout << "\nactual total interpolation point is" << n;
}

void generate_interpolation_zeros(dvector* r, int N, double length) {

    double dl, dl_2;

    double* zeros, _pi = 3.1415927;

    int n, i, j, k;

    zeros = new double[N];

    dl = length / (N - 1); // ÿ����ֵ�ڵ�֮��ļ��

    dl_2 = dl / 2; // �ڵ����һ��

    for (k = 0; k < N; k++)

        // zeros[k]=length/2.*(1-cos((2.*k+1.)/(2.*N)*_pi));  //chebyshev nodes

        // zeros[k]=length/2.*(1-(cos((2.*k+1.)/(2.*N)*_pi))/(cos(_pi/(2.*N)))); //Li's chebyshev grid

        // zeros[k]=length/2.*(1-cos(k*_pi/(N-1)));

        // zeros[k]=length/2.*(1+kesai[N-2][k]);//gauss integral point

        // zeros[k]=length/2.*(1+tan((k+0.5)/(2.*N)*_pi-_pi/4.)); //tangent nodes

        zeros[k] = k * dl; // uniform �����ֵ�ڵ�����꣬���� N �����ȷֲ��Ľڵ㣻�ڵ��ֵΪ k * dl���� 0 �� length ֮���� dl Ϊ���

    n = 0;

    for (i = 0; i < N; i++) {

        for (j = 0; j < N; j++) {

            for (k = 0; k < N; k++) {

                // r[n]=dvector(i*dl,j*dl,k*dl);

                r[n] = dvector(zeros[i], zeros[j], zeros[k]); // ������ά�Ľڵ�����

                n++;

            }
        }
    }

    for (k = 0; k < (N - 1); k++) {

        // zeros[k]=length/2.*(1-cos((2.*k+1.)/(2.*(N-1))*_pi)); //chebyshev nodes

        // zeros[k]=length/2.*(1-(cos((2.*k+1.)/(2.*(N-1))*_pi))/(cos(_pi/(2.*(N-1)))));//Li's chebyshev grid

        // zeros[k]=length/2.*(1+kesai[N-3][k]); //gauss integral point

        // zeros[k]=length/2.*(1+tan((k+0.5)/(2.*(N-1))*_pi-_pi/4.));    //tangent nodes

        zeros[k] = dl_2 + k * dl; // uniform
    }

    for (i = 0; i < (N - 1); i++) {
        for (j = 0; j < (N - 1); j++) {
            for (k = 0; k < (N - 1); k++) {
                // r[n]=dvector(dl_2+i*dl,dl_2+j*dl,dl_2+k*dl);
                r[n] = dvector(zeros[i], zeros[j], zeros[k]);
                n++;
            }
        }
    }
    cout << "Final n: " << n << endl;
    // ��ӡ��ά����
    std::cout << std::fixed << std::setprecision(4); // ���������ʽ
    std::cout << "Generated 3D coordinates:\n";
    for (int idx = 0; idx < n; idx++) {
        std::cout << "Point " << idx + 1 << ": (" << r[idx].vect_x << ", " << r[idx].vect_y << ", " << r[idx].vect_z
            << ")\n";
    }

    // ������̬������ڴ�
    delete[] zeros;
}

void generate_tartan_interpolation_zeros(dvector* r, int N, double length, double *zeroX, double *zeroY, double *zeroZ) {

    double dl, dl_2;

    double* zeros, _pi = 3.1415927;

    int n, i, j, k;

    zeros = new double[N];

    dl = length / (N - 1);

    dl_2 = dl / 2;

    for (k = 0; k < N; k++) {

        // zeros[k]=length/2.*(1-(cos((2.*k+1.)/(2.*N)*_pi))/(cos(_pi/(2.*N)))); //Li's chebyshev grid

        // zeros[k]=length/2.*(1-cos((2.*k+1.)/(2.*N)*_pi));//chebyshev node

        // zeros[k]=length/2.*(1-cos(k*_pi/(N-1)));

        // zeros[k]=length/2.*(1+kesai[N-2][k]);  //gauss integral point

        zeros[k] = k * dl; // uniform

        zeroX[k] = k * dl;
        zeroY[k] = k * dl;
        zeroZ[k] = k * dl;
    }

    n = 0;

    for (i = 0; i < N; i++)

        for (j = 0; j < N; j++)

            for (k = 0; k < N; k++) {

                // r[n]=dvector(i*dl,j*dl,k*dl);

                r[n] = dvector(zeros[i], zeros[j], zeros[k]);

                n++;
            }
}

void generate_GG(double _k, Complex* GG, dvector* r, dvector rs, int N) // only for staggered tartan grid

{
    int n, i, j, k;

    double R;

    n = 0;

    for (i = 0; i < N; i++)

        for (j = 0; j < N; j++)

            for (k = 0; k < N; k++) {

                R = (r[n] - rs).abs_vector();

                GG[n] = exp(Complex(0, 1) * _k * R) / R;

                n++;

            }

    for (i = 0; i < N - 1; i++)

        for (j = 0; j < N - 1; j++)

            for (k = 0; k < N - 1; k++) {

                R = (r[n] - rs).abs_vector();

                GG[n] = exp(Complex(0, 1) * _k * R) / R;

                n++;

            }

}

void generate_GG(int num_sp, double _k, Complex* GG, dvector* r, dvector rs) {

    int n;

    double R;

    for (n = 0; n < num_sp; n++) {

        R = (r[n] - rs).abs_vector();

        GG[n] = exp(Complex(0, 1) * _k * R) / R;
    }
}

void generate_the_Monomials_matrix_fortran_format(double* A, int N, double* spx, double* spy, double* spz, double c,
    double _k, double lambda) {

    int i, j;

    double* a, Ri;

    vector Ri_;

    for (i = 0; i < N; i++) {

        for (j = 0; j < N; j++) {

            Ri_ = vector((spx[i] - spx[j]), (spy[i] - spy[j]), (spz[i] - spz[j]));

            Ri = Ri_.abs_vector();

            //Ri = Ri / lambda; //+++++++++++++++++++++++++++++++++++++

            // if(Ri<1e-20)a[j]=0;

            // else

            // a[j]=Ri*Ri*log(Ri);

            // A[j*N+(i+1)]=1./sqrt(Ri*Ri+c*c);

            // A[j*N+(i+1)]=cos(_k*Ri)/sqrt(Ri*Ri+c*c);

            // A[j*N+(i+1)]=sqrt(Ri*Ri+c*c);

            // A[j*N+(i+1)]=exp(c*Ri*Ri);

            //     if(Ri<1e-20)

            //       A[i*N+(j+1)]=0;

            //		else

            //       A[i*N+(j+1)]=Ri*Ri*log(Ri);   //Thin plate spline

            //     if(Ri<1e-20)

            //       A[i*N+(j+1)]=0;

            //		else

            //       A[i*N+(j+1)]=Ri*Ri*Ri;   //Cubic

            //  A[i*N+(j+1)]=pow(sqrt(c*c+Ri*Ri),3);  //Natural cubic spline

            // A[i*N+(j+1)]=log(c*c+Ri*Ri); //Multilog

            //       A[i*N+(j+1)]=Ri;//linear

            // A[i*N+(j+1)]=cos(_k*Ri)/sqrt(Ri*Ri+c*c);

             A[i*N+(j+1)]=1./sqrt(Ri*Ri+c*c); //Inverse multiquadric

            // A[i*N+(j+1)]=exp(-c*Ri*Ri); //Gaussian

            // A[i*N+(j+1)]=1./(Ri*Ri+c*c);  //Inverse quadratic

            // A[i*N+(j+1)]=sqrt(Ri*Ri+c*c);//Multiquadric



            //if (Ri > 1) // with compact support

            //    A[i * N + (j + 1)] = 0;

            //else

            //    // A[i*N+(j+1)]=pow((1-Ri),4)*(4*Ri+1);    //C2

            //    A[i * N + (j + 1)] = pow((1 - Ri), 6) * (35 * Ri * Ri + 18 * Ri + 3); // C4
        }
    }
}

// R: the array used to store the interpolation points.

double Radial_Lagrange(dvector& r, dvector& ri, double* R, double R1234) {

    double f, Ri;

    dvector dr = r - ri;

    Ri = dr.abs_vector();

    // if(Ri>R[2])f=0;

    // else

    // f=1-Ri*Ri/R[1]*R[1];

    f = (R[1] - Ri) * (R[2] - Ri) / R[1] / R[2];

    return f;
}

double radial(dvector& r, dvector& ri, double c, double lambda) {

    double f;

    // Complex f;

    double Ri;

    dvector dr = r - ri;

    Ri = dr.abs_vector();

    Ri = Ri / lambda;

    // f=1./sqrt(Ri*Ri+c*c);  //Inverse multiquadric

    // f=exp(-c*Ri*Ri);  //Gaussian

    // f=1./(Ri*Ri+c*c);   //Inverse quadratic

    // f=sqrt(Ri*Ri+c*c);//Multiquadric

    // f=exp(Complex(0,1)*c*Ri);  //*Ri

    // f=log(Ri*Ri+c*c);  //shifted logarithm

    if (Ri > 1) // with compact support

        f = 0;

    else

        // f=pow((1-Ri),4)*(4*Ri+1);

        f = pow((1 - Ri), 6) * (35 * Ri * Ri + 18 * Ri + 3);

    // if(Ri<1e-20)f=0;

    // else

    // f=Ri*Ri*log(Ri);

    // f=pow(Ri,9.);

    // f=pow(10.,-c*Ri*Ri);

    // f=Ri*Ri*Ri;

    return f;
}

void scattered_point_lagrange_interpolation(Complex& Ga, double xx, double yy, double zz, Complex* GG, double* P[],
    double* a, int N, double c,
    double* spx, double* spy, double* spz, double _k, double lambda) {

    int n, p, q;

    double* t, temp;

    dvector Ri_;

    double Ri;

    for (n = 0; n < N; n++) {

        Ri_ = dvector(xx - spx[n], yy - spy[n], zz - spz[n]);

        Ri = Ri_.abs_vector();

        //Ri = Ri / lambda; //++++++++++++++++++++++++++++++++++++++++++++

        // if(Ri<1e-20)

        //  a[n]=0;

        // else

        //  a[n]=Ri*Ri*log(Ri);  //Thin plate spline

        // if(Ri<1e-20)

        //   a[n]=0;

        // else

        //	  a[n]=Ri*Ri*Ri;  //Cubic

        // a[n]=pow(sqrt(c*c+Ri*Ri),3);  //Natural cubic spline

        // a[n]=log(c*c+Ri*Ri); //Multilog

        //  a[n]=Ri;  //Linear

         a[n]=1/sqrt(Ri*Ri+c*c); //Inverse multiquadric

        // a[n]=cos(_k*Ri)/sqrt(Ri*Ri+c*c);

        // a[n]=sqrt(Ri*Ri+c*c); //Multiquadric

        // a[n]=exp(-c*Ri*Ri); //Gaussian

        // a[n]=1/(Ri*Ri+c*c);//Inverse quadratic



        //if (Ri > 1) // with compact support

        //    a[n] = 0;

        //else

        //    // a[n]=pow((1-Ri),4)*(4*Ri+1); //C2

        //    a[n] = pow((1 - Ri), 6) * (35 * Ri * Ri + 18 * Ri + 3); // C4
    }

    Ga = COMPLEX(0);

    for (p = 0; p < N; p++) {

        t = P[p];

        temp = 0;

        for (q = 0; q < N; q++) {

            temp = temp + t[q] * a[q];
        }

        Ga = Ga + temp * GG[p];
    }
}

/*void interpolation_test1()

{

   double lambda=1,cube_length=lambda/2.,distance;

   int num_p=3; //This is fixed.

   dvector *rr;

   int num_sca=num_p*num_p*num_p+(num_p-1)*(num_p-1)*(num_p-1);  //number of scattered points.

   rr=new dvector[num_sca];



   generate_interpolation_zeros(rr,num_p,cube_length);



   dvector rs; //The source point.

   distance=cube_length*2;

   rs=dvector(distance,cube_length/2.,cube_length/2.);



   dvector rf; //The field point;

   double xf,yf,zf;



   ddd:

   cout<<"\n input xf,yf, zf:";

   cin>>xf>>yf>>zf;

   rf=dvector(xf,yf,zf);



   double _k=2*3.1415927/lambda;



   COMPLEX G;

   double R=(rs-rf).abs_vector();

   G=exp(COMPLEX(0,1)*_k*R)/R;



   COMPLEX *GG;

   GG=new COMPLEX[num_sca];



   generate_GG(_k,GG,rr,rs,num_p);



   double RR[5],R1234;

   RR[1]=sqrt(3)*cube_length/4.;

   RR[2]=cube_length/2.;

   RR[3]=sqrt(2)*cube_length/2.;

   RR[4]=sqrt(3)*cube_length/2.;

   R1234=RR[1]*RR[2]*RR[3]*RR[4];



   COMPLEX f;

   double ff;

   int i;

   f=COMPLEX(0,0);

   for(i=0;i<num_sca;i++)

   {

      f=f+GG[i]*Radial_Lagrange(rf,rr[i],RR,R1234);

   }





   cout<<"\nf="<<f;

   cout<<"G="<<G;

   double error=abs((f-G)/G);

   cout<<"\nerror="<<error;

   goto ddd;



} */

//void interpolation_test2() {
//
//    double lambda = 1, cube_length = 4 * lambda, distance;
//
//    int num_p; // This is fixed.
//
//    dvector* rr;
//
//    int num_sca; // number of scattered points.
//
//    int num_sca0; // num_p*num_p*num_p;
//
//    ofstream fout;
//
//    fout.open("error_rad.txt", ios::ate);
//
//    if (!fout.is_open()) {
//        cerr << "Failed to open error_rad.txt" << endl;
//        return;
//    }
//
//    // generate_Gaussian_nodes();
//
//    cout << "\ninput 1D interpolation number: ";
//
//    cin >> num_p;
//
//    num_sca = num_p * num_p * num_p + (num_p - 1) * (num_p - 1) * (num_p - 1);
//
//    num_sca0 = num_sca;
//
//    rr = new dvector[num_sca];
//
//    Complex* GG;
//
//    GG = new Complex[num_sca];
//
//    matrixc16_ A, yy, gg;
//
//    A = matrixc16_(num_sca0, num_sca0);
//
//    yy = matrixc16_(num_sca0, 1);
//
//    gg = matrixc16_(num_sca0, 1);
//
//    //////////////////////////////////
//
//    /*   matrix_ A,yy,gg;
//
//       A=matrix_(num_sca0,num_sca0);
//
//       yy=matrix_(num_sca0,1);
//
//       gg=matrix_(num_sca0,1);*/
//
//       ///////////////////////////////////
//
//    generate_interpolation_zeros(rr, num_p, cube_length);
//
//    dvector rs; // The source point.
//
//    distance = 2 * cube_length;
//
//    rs = dvector(distance, cube_length / 2., cube_length / 2.);
//
//    dvector rf; // The field point;
//
//    double xf, yf, zf;
//
//    //   ddd:
//
//    //   cout<<"\n input xf,yf, zf:";
//
//    //   cin>>xf>>yf>>zf;
//
//    //   rf=dvector(xf,yf,zf);
//
//    double _k = 2 * 3.1415927 / lambda;
//
//    Complex G;
//
//    double R; //=(rs-rf).abs_vector();
//
//    //   G=exp(COMPLEX(0,1)*_k*R)/R;
//
//    // COMPLEX *GG;
//
//    // GG=new COMPLEX[num_sca];
//
//    generate_GG(_k, GG, rr, rs, num_p);
//
//    Complex f;
//
//    int i, j;
//
//    double c;
//
//    c = cube_length / (num_p - 1);
//
//    cout << "input parameter c: ";
//
//    cin >> c;
//
//    // c=sqrt(3)*c*1.5;
//
//    // c=6./c/c;
//
//    for (i = 0; i < num_sca0; i++)
//
//        gg[i] = GG[i];
//
//    // for()
//
//    //{
//
//    for (i = 0; i < num_sca0; i++)
//
//        for (j = 0; j < num_sca0; j++) {
//
//            A[i * num_sca0 + j] = Complex(radial(rr[i], rr[j], c, lambda));
//        }
//
//    // foutput(A,"A.txt");  //to check if A is in ill-condition
//
//    // cgmom(A,yy,gg,1e-6,0);
//
//    invert(A);
//
//    multi(yy, A, gg);
//
//    int I0, J0, K0, k, l;
//
//    I0 = 31;
//    J0 = 31;
//    K0 = 31;
//
//    double dll = cube_length / (I0 - 1);
//
//    double error = 0;
//
//    for (i = 0; i < I0; i++)
//
//        for (j = 0; j < J0; j++)
//
//            for (k = 0; k < K0; k++) {
//
//                rf = dvector(i * dll, j * dll, k * dll);
//
//                f = Complex(0);
//
//                for (l = 0; l < num_sca0; l++)
//
//                    f = f + yy[l] * Complex(radial(rf, rr[l], c, lambda));
//
//                R = (rs - rf).abs_vector();
//
//                G = exp(Complex(0, 1) * _k * R) / R;
//
//                error = error + abs((f - G) / G);
//            }
//
//    error /= I0 * J0 * K0;
//
//    cout << "error=" << error;
//
//    fout << num_sca0 << ' ' << "c=" << c << " error=" << error << '\n';
//
//    fout.close();
//
//    //   goto ddd;
//}

// ��ȡ��ǰʱ�䲢��ʽ��Ϊ�ַ���
std::string getCurrentTime() {
    // ��ȡ��ǰʱ��
    time_t now = time(nullptr);

    // ��ʽ��ʱ��
    struct tm tmBuffer;
    if (_localtime64_s(&tmBuffer, &now) != 0) { // ʹ�� _localtime64_s
        std::cerr << "Failed to get local time." << std::endl;
        return "";
    }

    char formattedTime[20];
    if (strftime(formattedTime, sizeof(formattedTime), "%Y-%m-%d %H:%M:%S", &tmBuffer) == 0) {
        std::cerr << "Failed to format time." << std::endl;
        return "";
    }

    return std::string(formattedTime);
}

// ����˷�����ʵ�֣�A��m*n��B��n*p��C��m*p
void matrixMultiply(const double* A, const double* B, double* C, int m, int n, int p) {
    // ��ʼ��������� C Ϊ 0
    for (int i = 0; i < m * p; ++i) {
        C[i] = 0.0;
    }

    // ���о���˷�
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < p; ++j) {
            for (int k = 0; k < n; ++k) {
                C[i * p + j] += A[i * n + k] * B[k * p + j];
            }
        }
    }
}

// ����ת�ú���ʵ�֣�A��m*n��B = A.t()
void matrixTranspose(const double* A, double* B, int m, int n) {
    // ���о���ת��
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            B[j * m + i] = A[i * n + j];
        }
    }
}

void new_test_scattered_lagrange(const double length_new, const double lambda, const double freq, ofstream& lxl_outfile, const string filename, const double c_new, std::deque<double> &Optimal_C_list);

void new_test_scattered_lagrange(const double length_new, const double lambda, const double freq, ofstream& lxl_outfile, const string filename, const double c_new, std::deque<double> &Optimal_C_list) {
     
    double pi = 3.1415926535897932384626433832795028841971693993751;

    int numsp_1D; // numsp_1D: the number of scattered points in 1 dimension.

    // double length=3,lambda=1;

    const double miu0 = 1.2566370614359172953850573533118e-6;
    const double epsilon0 = 8.85418781762039e-12;
    const double c_wave_speed = 1.0 / sqrt(miu0 * epsilon0);
    //double lambda = c_wave_speed / freq;
    //double lambda = 0.333333333;

    double length = length_new;

    //double lambda = 0.333333333;

    string output_path = filename;
    //ofstream lxl_outfile;
    lxl_outfile.open(output_path, std::ios_base::app);
    lxl_outfile << "X_length = " << length << "  Y_length = " << length << "  Z_length = " << length << "\n";
    lxl_outfile << "X_length(λ) = " << length / lambda << "  Y_length(λ) = " << length / lambda << "  Z_length(λ) = " << length / lambda << "\n";

    double c; // c: The shape parameter used in radial basis functions.

    int i, j, k;

    ///////////////////////////////////////

    // Generate the scattered points;

    ///////////////////////////////////////

    int test_type = 1;

    // test_type=1: Test the average interpolation error of the green's function generated

    // from the source point that is one cube length away from the field cube using the staggered tartan grid.

    // test_type=2: Test the average interpolation error of the green's function generated

    // from the source point that is one cube length away from the field cube using the ordinary tartan grid.

    // test_type=3: Test the average interpolation error of the green's function generated

    // from the source one cube length away from the field cube using the staggered tartan grid,

    // while the source point is moving from S to S' as in Fig. 7 of the paper.

    // test_type=4: Test the average interpolation error of the green's function generated

    // from the source one cube length away from the field cube using the ordinary tartan grid,

    // while the source point is moving from S to S' as in Fig. 7 of the paper.

    bool constant_c = true;

    // constant_c=false: the shape parameter c is varied.

    // constant_c=true: the shape parameter c is constant and should be input from the computer.

    cout << "\ntest_type=1: Test the average interpolation error of the green's function generated"

        << "\nfrom the source point that is one cube length away from the field cube using the staggered tartan grid."

        << "\ntest_type=2: Test the average interpolation error of the green's function generated"

        << "\nfrom the source point that is one cube length away from the field cube using the ordinary tartan grid."

        << "\ntest_type=3: Test the average interpolation error of the green's function generated"

        << "\nfrom the source one cube length away from the field cube using the staggered tartan grid,"

        << "\nwhile the source point is moving from S to S' as in Fig. 7 of the paper."

        << "\ntest_type=4: Test the average interpolation error of the green's function generated"

        << "\nfrom the source one cube length away from the field cube using the ordinary tartan grid,"

        << "\nwhile the source point is moving from S to S' as in Fig. 7 of the paper."

        << "\ntest_type=5: Test the average interpolation error of the green's function generated"

        << "\nfrom the source point that is one cube length away from the field cube using the new tartan grid.";

    //cout << "\n test_type=";

    //cin >> test_type;

    //cout << "\nplease input 1D interpolation number";

    //cin >> numsp_1D;

    //// cout<<"\nIs the parameter c constant(Yes:1|No:0)";

    //// cin>>constant_c;

    //cout << "\nplease input parameter c";

    //cin >> c;

    test_type = 2;
    numsp_1D = 3;
    c = 1.0;

    double IMQ_min_error = std::numeric_limits<double>::max();
    double IMQ_optimal_c = 0.0;
    Complex IMQ_optimal_G = Complex(0.0, 0.0);
    Complex Real_G = Complex(0.0, 0.0);
    bool Print_Interpolation_Point = false;
    double SVD_min_error = std::numeric_limits<double>::max();
    double SVD_max_error = std::numeric_limits<double>::min();

    c = c * length / lambda; //@@@@@@@@@@@@@@@

    // cout<<"\nThe number of intepolation point along one dimension is";

    // cin>>numsp_1D;

    int num_sp; // the total number of interpolation points of a tartan grid or a staggered one.

    if (test_type == 2 || test_type == 4)
        num_sp = numsp_1D * numsp_1D * numsp_1D;

    else if (test_type == 1 || test_type == 3)
        num_sp = numsp_1D * numsp_1D * numsp_1D + (numsp_1D - 1) * (numsp_1D - 1) * (numsp_1D - 1);

    else if (test_type == 5)
        num_sp = 8 * (numsp_1D - 1) * (numsp_1D - 1) * (numsp_1D - 1) + numsp_1D * numsp_1D * numsp_1D;

    cout << "\nnum_sp=" << num_sp;

    generate_Gaussian_nodes();

    // cout<<"\nThe wave length is(Note the cube length is unit):";

    // cin>>lambda;

    // length: the cube length.

    // lambda: the wavelength of the EM wave.

    dvector* rr; // Array used to store the interpolation points.

    double* spx, * spy, * spz;

    double* P[10000];

    double dc = 0.1;

    int nc = 0;

    ofstream fo;

    ofstream fot; // test file

    // fot.open("testfile.txt");

    Complex* GG;

    GG = new Complex[num_sp]; //*10

    double* SA;

    double* svdA, * svdV, * svdW; //,*P; //svdA=svdU*svdW*svdVT, after svd decomposition svdA is replaced by svdU

    SA = new double[num_sp * num_sp + 1];
    for (int i = 0; i < num_sp * num_sp + 1; i++) SA[i] = 0.0;

    svdA = new double[(num_sp + 1) * (num_sp + 1)];
    for (int i = 0; i < (num_sp + 1) * (num_sp + 1); i++) svdA[i] = 0.0;

    svdV = new double[(num_sp + 1) * (num_sp + 1)];
    for (int i = 0; i < (num_sp + 1) * (num_sp + 1); i++) svdV[i] = 0.0;

    svdW = new double[(num_sp + 1)];
    for (int i = 0; i < (num_sp + 1); i++) svdW[i] = 0.0;

    double* Recovered_A = new double[num_sp * num_sp];
    double* A_new = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++) {
        Recovered_A[i] = 0.0;
        A_new[i] = 0.0;
    }
    double* svd_U = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++) svd_U[i] = 0.0;
    double* svd_S = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++) svd_S[i] = 0.0;
    double* svd_V = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++) svd_V[i] = 0.0;

    for (i = 0; i < num_sp; i++)

        P[i] = new double[num_sp]; // Allocate memory for matrix A.

    rr = new dvector[num_sp]; // *10

    spx = new double[num_sp];

    spy = new double[num_sp];

    spz = new double[num_sp];

    double* zeroX = new double[numsp_1D];
    double* zeroY = new double[numsp_1D];
    double* zeroZ = new double[numsp_1D];

    double* poly;

    poly = new double[num_sp];

    double dlength;

    dlength = 0.1;

    int NC;

    if (test_type == 1 || test_type == 2 || test_type == 5)
        NC = 20;

    else
        NC = 15;


    double max_c = 0.0;
    double c_step = 0.1;
    if (c_new < 0.0) {
        max_c = 10.0;
        c = 0.1;
    }
    else {
        max_c = c_new;
        c = c_new;
    }
    //for (int nc = 0; nc < NC; nc++) {
    for(; c <= max_c; c += c_step) {
            // cout<<"\nnc="<<nc;

            // nc=30;

            // if(test_type==3||test_type==4||constant_c)

            //{

            //	cout<<"\n Please input a constant shape parameter C:";

            //   cin>>c;

            //}

            // else

            // c=(0.5+nc*dc)*length/(numsp_1D-1);  //c: used in radial basis functions.

            // c=length/(numsp_1D-1);

            //c = 0.1 + 0.1 * nc;

        if (test_type == 1 || test_type == 3)

            generate_interpolation_zeros(rr, numsp_1D, length);

        else if (test_type == 2 || test_type == 4)

            generate_tartan_interpolation_zeros(rr, numsp_1D, length, zeroX, zeroY, zeroZ);

        else if (test_type == 5)

            generate_new_interpolation_zeros(rr, numsp_1D, length);

        for (i = 0; i < num_sp; i++) {

            spx[i] = rr[i].vect_x;

            spy[i] = rr[i].vect_y;

            spz[i] = rr[i].vect_z;
        }
           
        /*if (!Print_Interpolation_Point) {
            for (i = 0; i < num_sp; i++) {
                lxl_outfile << "x=" << spx[i] << "  y=" << spy[i] << "  z=" << spz[i] << "\n";
            }
            Print_Interpolation_Point = true;
        }*/
        

        //////////////////////////////////////////////////////////////

        // Generate the Green's functions and store them in array GG.

        //////////////////////////////////////////////////////////////

        double _k = 2 * pi / lambda;

        double distance = 2 * length;

        // distance:The distance from the original point to the nearest point of the nearest far group along x axis.

        dvector rs; // Source point.

        if (test_type == 1 || test_type == 2 || test_type == 5)

            rs = shift + dvector(distance, 0.5 * length, 0.5 * length);

        else
            rs = shift + dvector(distance, (0.5) * length, (0.5 + (dlength * nc)) * length);

        //lxl_outfile << "\nrs = (" << rs.vect_x << ", " << rs.vect_y << ", " << rs.vect_z << ")\n";

        generate_GG(num_sp, _k, GG, rr, rs);//rr�洢�����еĲ�ֵ�㣬�������������rs�����еĲ�ֵ��֮������и��ֺ���Ȼ�󱣴�����GG��

        // generate_GG(_k,GG,rr,rs,numsp_1D);

        // for(i=0;i<num_sp;i++)

        //{

        //   fot<<i<<" Green="<<real(GG[i])<<endl;

        //}

        ////////////////////////////////////////

        // Generate the Monomials matrix.

        ////////////////////////////////////////

        // for(i=0;i<num_sp;i++)

        //{

        //	A[i]=new double[num_sp];  //Allocate memory for matrix A.

        //}

        // generate_the_Monomials_matrix(A,num_sp,spx,spy,spz);

        // generate_the_Monomials_matrix(A,num_sp,spx,spy,spz,c);

        // generate_the_Monomials_matrix(AA,num_sp,spx,spy,spz,c);

        // double *GINVA=new double[num_sp*num_sp+1],

        // double TOL=1e-10;

        // int IRANK;


        generate_the_Monomials_matrix_fortran_format(SA, num_sp, spx, spy, spz, c, _k, lambda);

        /*lxl_outfile << "\nSA:" << "\n";
        int SA_size = num_sp * num_sp + 1;
        for (int i = 0; i < SA_size; i++) {
            lxl_outfile << std::setprecision(16) << SA[i] << "  ";
            if (i % num_sp == 0) lxl_outfile << "\n";
        }*/
       

        /*   k=0;

           for(i=0;i<num_sp;i++)

           {         for(j=0;j<num_sp;j++)

                    {

                      fot<<"A"<<k<<"="<<SA[i*num_sp+(j+1)]<<"\t";

                      k++;

                    }

                    fot<<endl;

           } */

        for (i = 0; i < num_sp; i++)

            for (j = 0; j < num_sp; j++) {

                svdA[(i + 1) * (num_sp + 1) + j + 1] = SA[j * num_sp + (i + 1)];
            }

        //lxl_outfile << "\n\nsvdA:" << "\n";
        for (int i = 0; i < num_sp; i++) {
            //lxl_outfile << "{";
            for (int j = 0; j < num_sp; j++) {             
                A_new[i * num_sp + j] = svdA[(i + 1) * (num_sp + 1) + j + 1];
                //lxl_outfile << std::setprecision(16) << svdA[(i + 1) * (num_sp + 1) + j + 1];
                //if (j < num_sp - 1) lxl_outfile << ",  ";
            }
            //lxl_outfile << "},\n";
        }

        svdcmp(svdA, num_sp, num_sp, svdW, svdV);


        //lxl_outfile << "\n\nsvd_U:" << "\n";
        for (int i = 0; i < num_sp; i++) {
            for (int j = 0; j < num_sp; j++) {
                svd_U[i * num_sp + j] = svdA[(i + 1) * (num_sp + 1) + j + 1];
                //lxl_outfile << std::setprecision(16) << svdA[(i + 1) * (num_sp + 1) + j + 1] << "   ";
            }
            //lxl_outfile << "\n";
        }
        //lxl_outfile << "\n";


        //lxl_outfile << "\n\nsvd_S:" << "\n";
        for (int i = 0; i < num_sp; i++) {
            svd_S[i * num_sp + i] = svdW[i + 1];
            //lxl_outfile << std::setprecision(16) << svdW[i + 1] << "   ";
        }
        //lxl_outfile << "\n";


        //lxl_outfile << "\n\nsvd_V:" << "\n";
        for (int i = 0; i < num_sp; i++) {
            for (int j = 0; j < num_sp; j++) {
                svd_V[i * num_sp + j] = svdV[(i + 1) * (num_sp + 1) + j + 1];
                //lxl_outfile << std::setprecision(16) << svdV[(i + 1) * (num_sp + 1) + j + 1] << "   ";
            }
            //lxl_outfile << "\n";
        }
        //lxl_outfile << "\n";

        //ʵ��A����ĸ�ԭ A = U * S * V.t()
        int m = num_sp, n = num_sp, p = num_sp;
        //lxl_outfile << "Recover A = U * S * V.t()";
        double* US = new double[m * n];
        matrixMultiply(svd_U, svd_S, US, m, m, n);
        double* svd_Vt = new double[n * n];
        matrixTranspose(svd_V, svd_Vt, n, n);
        matrixMultiply(US, svd_Vt, Recovered_A, m, n, n);


       /* lxl_outfile << "\nsvd_Vt:" << "\n";
        for (int i = 0; i < num_sp * num_sp; i++) {
            lxl_outfile << std::setprecision(16) << svd_Vt[i] << "   ";
            if ((i + 1) % num_sp == 0) {
                lxl_outfile << "\n";
            }
        }
        lxl_outfile << "\n";*/


        /*lxl_outfile << "\n\nRecovered_A : \n";
        for (int i = 0; i < num_sp * num_sp; i++) {
            lxl_outfile << std::setprecision(16) << Recovered_A[i] << "   ";
            if ((i + 1) % num_sp == 0) lxl_outfile << "\n";
        }
        lxl_outfile << "\n";*/
          
        double Recover_Avergae_error = 0.0;
        for (int i = 0; i < num_sp * num_sp; i++) {
            Recover_Avergae_error += abs(A_new[i] - Recovered_A[i]) / abs(A_new[i]);
        }
        Recover_Avergae_error = Recover_Avergae_error / (num_sp * num_sp);
        if ((SVD_min_error - Recover_Avergae_error) > 1E-20) {
            SVD_min_error = Recover_Avergae_error;
        }
        if ((Recover_Avergae_error - SVD_max_error) > 1E-20) {
            SVD_max_error = Recover_Avergae_error;
        }

        ////////////////////////////////////////////////

        // Generate the polynomial coefficients.

        ////////////////////////////////////////////////

        //   SVD(&num_sp,&num_sp, A[1],&num_sp,&TOL,&IRANK,&GINVA[1],&num_sp);

        //   for(i=0;i<num_sp;i++)  //Initialize the value to the polynomials.

        //   for(j=0;j<num_sp;j++)

        //   {

        //   	P[i][j]=GINVA[(i+1)+num_sp*j];

        //   }

        for (i = 0; i < num_sp; i++)

            for (j = 0; j < num_sp; j++) {

                P[i][j] = 0;

                for (k = 0; k < num_sp; k++)

                    P[i][j] += svdV[(i + 1) * (num_sp + 1) + (k + 1)] * svdA[(j + 1) * (num_sp + 1) + (k + 1)] /
                    svdW[k + 1];
            }

        /*lxl_outfile << "\n\nP:" << "\n";
        for (int i = 0; i < num_sp; i++) {
            for (int j = 0; j < num_sp; j++) {
                lxl_outfile << std::setprecision(16) << P[i][j] << "   ";
            }
            lxl_outfile << "\n";
        }
        lxl_outfile << "\n";*/
        //lxl_outfile << "Recover_Avergae_error = " << std::setprecision(16) << Recover_Avergae_error << "\n";


        ///////////////////////////////////////////////////////////////////////////////

        // Test the accuracy of the interpolation in some uniformly spaced points.

        ///////////////////////////////////////////////////////////////////////////////

        int num_p = 4; // numsp_1D*2*2; //num_p: the number of test points in one dimension.

        double dl = length / (num_p - 1);

        double xx, yy, zz; // Used to store the coordinates of the current points.

        Complex G0, Ga; // G0: the true value; Ga: the approximated value.

        dvector rf;

        double R, error, errorsum;

        double errormin = 1e10, errormax = -1e10;

        errorsum = 0;

        double Test_x = length;
        double Test_y = zeroY[0] + (zeroY[1] - zeroY[0]) / 2.0;
        double Test_z = zeroZ[0] + (zeroZ[1] - zeroZ[0]) / 2.0;

        dvector rf_new = dvector(Test_x, Test_y, Test_z);
        //lxl_outfile << "rf_new = (" << rf_new.vect_x << ", " << rf_new.vect_y << ", " << rf_new.vect_z << ")\n";

        double R_new = (rs - rf_new).abs_vector();
        Real_G = exp(Complex(0, 1) * _k * R_new) / R_new;
        Complex RBF_G = Complex(0.0, 0.0);
        scattered_point_lagrange_interpolation(RBF_G, Test_x, Test_y, Test_z, GG, P, poly, num_sp, c, spx, spy, spz, _k, lambda);

        double error_new = abs(Real_G - RBF_G) / abs(Real_G);
        /*lxl_outfile << "c = " << std::setprecision(16) << c << "\n";
        lxl_outfile << "error = " << error_new << "\n";
        lxl_outfile << "Real_G = " << Real_G << "   RBF_G = " << RBF_G << "\n";*/
        
        if (IMQ_min_error - error_new > 1E-10) {
            IMQ_min_error = error_new;
            IMQ_optimal_c = c;
            IMQ_optimal_G = RBF_G;
        }

        //for (i = 0; i < num_p; i++)

        //    for (j = 0; j < num_p; j++)

        //        for (k = 0; k < num_p; k++) {

        //            xx = i * dl;
        //            yy = j * dl;
        //            zz = k * dl;
        //            lxl_outfile << "(xx, yy, zz) = " << "(" << xx << ", " << yy << ", " << zz << ",)\n";

        //            rf = dvector(xx, yy, zz);

        //            R = (rs - rf).abs_vector();

        //            G0 = exp(Complex(0, 1) * _k * R) / R;

        //            scattered_point_lagrange_interpolation(Ga, xx, yy, zz, GG, P, poly, num_sp, c, spx, spy, spz, _k,
        //                lambda);

        //            // scattered_point_lagrange_interpolation(Ga,xx,yy,zz,GG,P,poly,num_sp);

        //            error = abs((G0 - Ga) / G0);
        //            lxl_outfile << "error = " << error << "   G0 = " << G0 << "   Ga = " << Ga << "\n";
        //            lxl_outfile << "----------------------------------------------------\n\n";

        //            if (errormin > error)
        //                errormin = error;

        //            if (errormax < error)
        //                errormax = error;

        //            // cout<<"\n error ="<<error;

        //            errorsum += error;
        //        }

        //error = errorsum / (num_p * num_p * num_p);

        //lxl_outfile << "\n Minimum error =" << errormin << "\n";

        //lxl_outfile << "\n Maximum error =" << errormax << "\n";

        //lxl_outfile << "\n Average error =" << error << "\n";

        //lxl_outfile << "----------------------------------------------------------------\n";

       /* fo.open("error_rad_svd.txt", ios::ate);

        if (error) {

            if (test_type == 3 || test_type == 4)

                fo << num_sp << " d=" << length / lambda << " c=" << c << " souce position="
                << (0.5 + nc * dlength) * length << " error=" << error << '\n';

            else

                fo << num_sp << " d=" << length / lambda << " c= " << c << " error= " << error << '\n';
        }

        fo.close();*/

        //   if((test_type==1||test_type==2)&&constant_c)

        //   {

        //   	break;

        //   }
    }
    Optimal_C_list.push_back(IMQ_optimal_c);

    lxl_outfile << "Real_G = " << std::setprecision(18) << Real_G << "\n";
    lxl_outfile << "IMQ_min_error = " << std::setprecision(18) << IMQ_min_error << "   ";
    lxl_outfile << "IMQ_optimal_c = " << std::setprecision(18) << IMQ_optimal_c << "   ";
    lxl_outfile << "IMQ_optimal_G = " << std::setprecision(18) << IMQ_optimal_G << "\n";
    //lxl_outfile << "\nsvd_min_error = " << std::setprecision(18) << SVD_min_error << "   svd_max_error = " << std::setprecision(18) << SVD_max_error << "\n";
    lxl_outfile << "=============================================================\n\n";

    lxl_outfile.close();
    delete[] zeroX;
    delete[] zeroY;
    delete[] zeroZ;
}

int main() {

    //    interpolation_test2(); // use invert function to calculate polynomial coefficients

    
    string output_path = "lxl_test_Reproduced_WHG.log";
    ofstream lxl_outfile(output_path, std::ios_base::trunc);
    string formattedTime = getCurrentTime();
    lxl_outfile << formattedTime << "\n";
    lxl_outfile.close();

    const double miu0 = 1.2566370614359172953850573533118e-6;
    const double epsilon0 = 8.85418781762039e-12;
    const double c_wave_speed = 1.0 / sqrt(miu0 * epsilon0);
    double freq = 15E+09;
    const double lambda = c_wave_speed / freq;
    //const double lambda = 0.333333333;
    double length = 0.25 * lambda;
    std::deque<double> Optimal_C_list, Optimal_C_list_useless;

    lxl_outfile.open(output_path, std::ios_base::app);
    lxl_outfile << "///  �̶� c == 0.55���� == 1.0���ֱ��� 2�� ~ 0.1�� ��Cube, step = 0.1�� ///\n";
    lxl_outfile.close();
    for (double i = 2.0; i > 0.05; i-=0.1) {
        double length = i;
        new_test_scattered_lagrange(length, 1.0, freq, lxl_outfile, output_path, 0.550, Optimal_C_list_useless);
    }

    lxl_outfile.open(output_path, std::ios_base::app);
    lxl_outfile << "///  �̶� c == 0.55���� == 1.0���ֱ��� 0.1�� ~ 0.01�� ��Cube, step = 0.01�� ///\n";
    lxl_outfile.close();
    for (double i = 0.10; i > 0.005; i -= 0.01) {
        double length = i;
        new_test_scattered_lagrange(length, 1.0, freq, lxl_outfile, output_path, 0.550, Optimal_C_list_useless);
    }

    exit(0);


    for (int i = 0; i < 8; i++) {
        length = (0.25 + i * 0.25) * lambda;
        new_test_scattered_lagrange(length, lambda, freq, lxl_outfile, output_path, -1, Optimal_C_list);
    }

    lxl_outfile.open(output_path, std::ios_base::app);
    lxl_outfile << "///  ʹ��2�˵�cube������cֵ���μ���1.95�ˡ�1.9�ˡ�1.85�ˡ�1.8�˵�cube ///\n";
    lxl_outfile.close();

    double c_new = 0.0;
    c_new = Optimal_C_list[7];
    c_new = 1.0;
    new_test_scattered_lagrange(2.0 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(1.95 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(1.9 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(1.85 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(1.8 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);

    lxl_outfile.open(output_path, std::ios_base::app);
    lxl_outfile << "///  ʹ��1.5�˵�cube������cֵ���μ���1.45�ˡ�1.4�ˡ�1.35�ˡ�1.3�˵�cube ///\n";
    lxl_outfile.close();

    c_new = Optimal_C_list[5];
    c_new = 1.5;
    new_test_scattered_lagrange(1.5 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(1.45 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(1.4 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(1.35 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(1.3 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);

    lxl_outfile.open(output_path, std::ios_base::app);
    lxl_outfile << "///  ʹ��1.0�˵�cube������cֵ���μ���0.99�ˡ�0.98�ˡ�0.97�ˡ�0.96�ˡ�0.95�˵�cube ///\n";
    lxl_outfile.close();

    c_new = Optimal_C_list[3];
    c_new = 2.0;
    new_test_scattered_lagrange(1.0 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(0.99 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(0.98 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(0.97 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(0.96 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
    new_test_scattered_lagrange(0.95 * lambda, lambda, freq, lxl_outfile, output_path, c_new, Optimal_C_list_useless);
 
    std::cout << "\nfinished!" << std::endl;

    //new_test_scattered_lagrange(); //use SVD function to calculate polynomial coefficients
}

#endif
